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We study theoretically the current voltage characteristics of intrinsic Josephson junctions in high- 
T c superconductors. An oscillation of the breakpoint current on the outermost branch as a function 
of coupling a and dissipation /3 parameters is found. We explain this oscillation as a result of the cre- 
ation of longitudinal plasma waves at the breakpoint with different wave numbers. We demonstrate 
the commensurability effect and predict a group behavior of the current-voltage characteristics for 
the stacks with a different number of junctions. A method to determine the wave number of longi- 
tudinal plasma waves from a- and /9-dependence of the breakpoint current is suggested. We model 
the a- and /9-dependence of the breakpoint current and obtain good agreement with the results of 
simulation. 



Creating new materials with given properties is an ac- 
tual problem of physics, chemistry, and material science. 
This is related to the system of Josephson junctions, too, 
which is a perspective object for superconducting elec- 
tronics and is being investigated intensively now. A sim- 
ulation of the current-voltage characteristics (IVC) of a 
stacks of intrinsic Josephson junctions (IJJ)i at different 
values of the model parameters such as the coupling and 
dissipation parameters is a way to predict the properties 
of the IJJ. McCumber and Steward have investigated the 
return current as a function of dissipation parameter in 
a single Josephson junction a long time ago^ In the case 
of the system of junctions, the situation is cardinally dif- 
ferent. The IVC of IJJ is characterized by a multiple 
branch structure and branches have a breakpoint region 
with its breakpoint current (BPC) and transition current 
to another branch^ The BPC is determined by the cre- 
ation of the longitudinal plasma waves (LPW) with a 
definite wave number k, which depends on the parame- 
ters a and /?, the number of junctions in the stack, and 
boundary conditions. If we neglect the coupling between 
junctions, the branch structure disappears, and the BPC 
coincides with the return current. As we know, an inves- 
tigation of the McCumber-Steward dependence for the 
different branches of IVC for IJJ has not been done yet. 
Machida and Koyama£ have stressed that capacitive cou- 
pling takes various values in HTSC and layered organic 
superconductors and they presented a systematic study 
for the capacitively coupled Josephson junctions (CCJJ) 
model, focusing on the dependence of phase dynamics 
on the strength of the capacitive coupling constant from 
weak to strong coupling regimes. But they did not in- 
vestgate the breakpoint region in the simulated IVC. 

In this Letter, we generalized the McCumber-Steward 
dependence of the return current for the case of IJJ in 
the HTSC. We investigate the BPC Ib P on the outermost 
branch as a function of the coupling a and dissipation 
parameters for the stacks with a different number of IJJ 
and demonstrate a plateau with BPC oscillation. Based 
on the idea of the parametric resonance in the stack of 
IJJ, a modeling of the cv/3-dependence of the BPC has 



been done, and good qualitative agreement with the re- 
sults of simulation has been obtained. We show that the 
a/3-dcpendcncc of the BPC is an instrument to determine 
the mode of LPW created at the breakpoint in the stacks 
with a different number of junctions. 

A system of dynamical equations in the capacitively 
coupled Josephson junctions model with diffusion current 
(CCJJ+DC model)&2 

d 2 , dipi . , . 

— ipi = (i - sin tpi - P—jj-) + a(sm <p l+1 + sin tp l _ 1 

-2mn<p l )+aft— + — -2—) (1) 

for the gauge-invariant phase differences (fi(t) = 0;+i(t) — 

&i(t) — x Ji +1 dzA z (z, t) between superconducting layers 
(5-layers) for the stacks with a different number of in- 
trinsic junctions has been numerically solved. Here 9i is 
the phase of the order parameter in S-layer I, A z is the 
vector potential in the barrier. 

The CCJJ+DC model is different from the CCJJ 
model&&i£ by the last term on the right hand side. This 
coupled Ohmic dissipation term might be derived by the 
microscopic theory^ or phenomenologically by the inclu- 
sion of the diffusion current between S-layers and leads to 
the equidistant branch structure in the IVCi 7 - The details 
concerning the system ([T]) are presented in Ref. 7 Here 
we use the periodic boundary conditions considering the 
first S-layer as a neighbor of the last one. 

The simulated IVC have the breakpoint on their out- 
ermost branches. We have calculated the /3-dependcncc 
of the BPC hp at fixed value of a, changing (3 in the 
interval (0,1) by step 0.005. The result of the calculation 
at a = 0, 1 and 5 is presented in Fig. QJ,. At a = 0, the 
IVC does not manifest the multibranch structure, and the 
breakpoint coincides with the return current. The curves 
at a 7^ have new features in comparison with the case 
without coupling. Particularly, they show a stronger in- 
crease of the Ib p at small (3, a plateau at Ib p — 0.83 and 
the oscillation of the Ib p on this plateau, and a transition 
to the non-hysteretic regime ( second plateau) at smaller 
[3. These features are discussed below. We change the 
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FIG. 1: (Color online) (a) - The /3-dependence of the BPC 
It p of the outermost branch in the IVC at different values of 
coupling parameter a; b) - The a/3-dependence of the for 
a stack of 10 IJJ. 



coupling parameter a in the interval (0,8) by step 0.1 
and repeat the calculations of the /3-dependence of hp. 
By this method, we build the three-dimensional picture 
of the a/3-dependence of the Ib p for a stack with 10 IJJ, 
which is shown in Fig. [Tp. We see two plateaus on this 
dependence and the oscillations of the Ib p on the first 
one as a function of a and /3. We note the next features 
for the /3-dependence : i) At a equal to zero, our results 
for /3-dependence of the h p coincide with the previous 
simulation of the /3-dependence of the return current 2 ; 
ii) at small /?, the /3-dependence is getting sharper with 
the increase in a; iii) the oscillations of the Ib p are get- 
ting stronger at larger a; iiii) with the increase in a, 
the transition to the non-hysteretic regime (to the sec- 
ond plateau) is approached at smaller /3 . For the in- 
dependence of the Ib p we may note: i) At small /3, the 
a-dependence is monotonic, and Ib p is increasing with a; 
ii) at some /3, the oscillations of Ib p appear, iii) with the 
increase in /3, the transition to the non-hysteretic regime 
is observed at smaller a. The value of the Ib p changes 
strongly at small a and [3. On the first plateau, the vari- 
ation of the Ib p consists of ~ 3 ~ 4 percent of the value 
of I c for N = 10. As we can see below, it depends on the 
number of junctions in the stack and decreases with N. 

Let us analyze in more detail the a- and /3-dependence 
of the Ib p ■ Fig. [TJi demonstrates the general features of /3- 
dependence of the Ib P at different values of the coupling 
parameter. To clearly show these features, we demon- 
strate in Fig. [2^ in an increased scale the /3-dependence 
of the Ib p at a = 3. We can see clearly four maxi- 
mums of Ib p on this curve. Using the Maxwell equation 
div(E/d) = 47r j o, we express the charge pi on the super- 
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FIG. 2: (Color online) a) - The /3-dependence of the Ib p for 
a stack with 10 IJJ at a — 3; b) - The a-dependence of the 
Ibp at j3 = 0.3; c) - The charge distribution among the layers 
corresponding to the different plasma modes in the stack of 
10 IJJ at a = 3 and p = 0.24, 0.27, 0.3, 0.4. 



conducting layer i by the voltages and V^j+i in 

the neighbor insulating layers p l = 4^3(^,1+1 - 
Solution of the system of equations (p} gives us the volt- 
ages V^j+i in all junctions in the stack, and it allows 
us to investigate the time dependence of the charge on 
each S-layer. We analyze the time dependence of the 
charge oscillations on S-layers at (3 equal to 0.24, 0.27, 
0.3 and 0.4 (around each maxima). The charge distri- 
butions among the S-layers in the stack at a fixed time 
moment at the breakpoint of the outermost branch are 
presented in Fig. [2J:. The charge oscillations on S-layers 
correspond to standing LPW with k equal to tt, 4tt/5, 
37r/5 and 2tt/5, relating to the four different intervals 
of the /3 with four maximums in this region. Fig. [2h 
shows the a-dependence of Ib p at /3 = 0.3, and it demon- 
strates four regions corresponding to the different modes 
of LPW. 

To prove our results and test the idea that at the break- 
point a parametric resonance is approached and plasma 
mode is excited by Josephson oscillations, we have mod- 
eled the a/3-dependence of the Ib p in the CCJJ+DC 
model. The equation for the Fourier component of the 
difference of phase differences 5(pi = tfi+ij — ¥1,1-1 be- 
tween neighbor junctions is 3 Sk+P{k)6k+cos(D,(k)T)Sk = 
0, where t = uj p (k)t, ui p (k) — tu p C, /3(fc) = j3C, 
fi(fc) = Sl/C and C = y/l + 2a(l - cos(fc)). This equa- 
tion shows a resonance with changing its parameters 
/3(A)) and £l(k). In Fig. [3^,, we have plotted the para- 
metric resonance region for this equation on the dia- 
gram /3(fc) — f2(fc). Using this diagram, we determine 
the curve which corresponds to the edge of the reso- 
nance region. This curve is shown in Fig. [3^i by dots. 
We consider that the point on this curve correspond- 
ing to maxQ(k) at a fixed value of /3(fc) gives us the 
value of the ttb P (k) which corresponds to the breakpoint 
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FIG. 3: (Color online) a) - Parametric resonance region in 
Q(k)~ (3(k) diagram. The value £l(k) — £lb P (k) corresponds to 
the breakpoint voltage on the outermost branch; b) - Result of 
modeling of the a/3-dependence of the h P for plasma modes 
with k = 7T and k = 27r/5 for a stack of 10 IJJ; c) - The 
modeled a-dependence of Ib p for stack with 10 IJJ at /? = 0.3 
corresponding to the creation of the LPW with different k; d) 
- The modeled /3-dependence of Ib p at a = 3. 



voltage. Taking into account the relations for the out- 
ermost branch flbp{k) = Vb p / (Ny^l + 2a(l — cos k)) and 
V bp /N = hp/ (3, we get 



/bp (a, 0, k) = Py/l + 2a(l - coB*)n^(fc,/3). (2) 

As an example, using the expression ^ for /ft p , we have 
plotted in FigOj the three-dimensional a/?-dependence of 
the Ibp for two plasma modes with k = ir and k — 27r/5 
for a stack with 10 IJJ. Comparing Fig[3)D with Fig[T|3, 
we note that the main features of the simulated and mod- 
eled a/3-dependence of the Ibp are in agreement. Using 
the formulas ([2|) , we have calculated the a-dependence of 
the Ibp at j3 = 0.3 for plasma modes with different wave 
numbers k. The corresponding curves are presented in 
Fig. [3J:. We see that these results of modeling coincide 
as well qualitatively with the results of simulation pre- 
sented Fig[T|D. Both kinds of curves show the same be- 
havior. We can see the increase in the distance between 
the maximums of Ibp and their sloping with increase in 
k in simulated and modeled curves. Fig. shows the 
modeled /3-dependence of Ibp at a = 3, which is obtained 
from the resonance region data. This dependence is in 
agreement with the results of simulation as well, and it 
demonstrates the oscillations of the Ib p , but it does not 
reflect the decrease in the values of Ibp maximums which 
is shown in Fig. [5^. This is a result of the approxima- 
tions we have used to obtain the linearized equation for 
the Fourier component of the difference of phase differ- 
ences for neighbor junctions^. The theoretical considera- 
tions which we use to model the a(3- dependence of the 



Ibp lead to the conclusion that there are regions on the 
a/3-dependence of Ibp which correspond to the creation 
of the LPW with a different wave number k and explain 
the origin of the h p oscillations. The ideas and results 
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FIG. 4: (Color online) The simulated IVC of the outermost 
branch in the stacks with a different number of junctions at 
a = 3, 13 = 0.3. 
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FIG. 5: (Color online) a) - The simulated /3-dependence of the 
h P for the stacks with 3, 6, 9 and 12 IJJ at a — 3. The region 
corresponding to the creation of the LPW mode with wave 
number k = 5tv/6 is shown by arrows . b) - The simulated 
a-dependence of the Ib P for the stacks with 5, 10 and 15 IJJ 
at p = 0.3. 

presented above have strong support from the results of 
investigation of the a- and /3-dependence of the Ibp in 
the case of a different number of IJJ in the stack. The 
minimal wavelength A which might be realized in the 
discrete lattice at periodic boundary conditions is two 
lattice units. So, in the stack with N junctions, the LPW 
with k = 27m /iV may exist, where n is an integer from 
1 to N/2 for even N and from 1 to (TV - l)/2 for odd 
N. Because of the term (1 — cosfc) in the LPW 

with k corresponding to the highest Ibp in the decreasing 
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current process is created. In RefQ], we showed that, at 
small values of a and (3 at periodic boundary conditions 
for stacks with even N, the 7r-mode of LPW is created, 
but for stacks with odd N the LPW with k = (N-l)n/N 
is observed. Here we consider a case of strong coupling 
between junctions, and the results are different from the 
previous consideration. 

Fig. Q] shows the result of simulation of the outermost 
branch in the IVC near the breakpoint for a stack with 
a = 3, j3 = 0.3 and N from N = 3 to N = 15. We can see 
that the value of Ib p depends on the number N of IJ J in 
the stack, excluding the stack with N = 3n, where n is an 
integer number. Time dependence analysis of the charge 
oscillations on the S-layers shows that, at the breakpoint 
in the stacks with N — 3n, the LPW with k = 2tt/3 is 
created. In the stack with N = 4, we observe the LPW 
with A = 4. We will not touch the question concerning 
the breakpoint region in the IVC presented in Fig. [U It 
will be considered in detail somewhere else. We may note 
another interesting group behavior of the IVC, presented 
in Fig. |U There is a monotonic increase of the It p with N 
for stacks with N = 3n + 1, n > 1. The same monotonic 
behavior was observed for stacks with N — 3n+2. Below, 
we explain these results using the idea of LPW creation 
at the breakpoint. 

Comparison of the a- or /3-dependence of the Ib p for 
stacks with a different number of IJJ give us a simple 
method to determine the wave numbers k of the LPW. 
Fig. shows the /3-dependence of the hp at a — 3 for 
the stacks with 3,6,9 and 12 IJJ. It demonstrates that, in 
some intervals of /3, the stacks with different N have the 
equal value of the Ib p - Particularly, all stacks have the 
equal values of the Ib p in some interval around (3 = 0.3. 
According to the results of modeling, for the stack with 
given N, the intervals on the the curves of the a- and 
/3-dependence, corresponding to the different modes of 
the LPW, follow in decreasing order in k. Because this 
interval around (3 — 0.3 corresponds to the region around 
the maximum on the /3-dependence of the Ibp for stack 
with N = 3, the second maximum for the stacks with 
N = 6 and N — 9, and the third maximum for the 
stack with N = 12, we may conclude that in this inter- 
val the LPW with k — 2ir/3 is created. For stacks with 
N = 6 this interval is continued until (3 = 0.365. Using 
this method of the wave number determination, which we 
call k — a/3-method, we can determine all modes of LPW 
which might be created in stacks with different parame- 
ters a and (3 and a different number of IJJ. Particularly, 



we find that, on the /3-dependence, the interval (0, 0.27) 
and the region (3 > 0.41 correspond to the creation of 
the 7T- and tt/3- modes of LPW, respectively. From the 
a-dependence of the Ib p which is presented in Fig. [SJa for 
stacks with 5, 10 and 15 IJJ, we find that the interval 
around the maximum with 2.35 and the region a > 4.82 
correspond to the creation of the 47r/5- and 7r/5- modes 
of LPW, respectively. 

Using the k — a/3-method, we find the values of k for 
IVC presented in Fig. 2J In the stacks with N = 3n 
(dash-dotted curves in Fig. Ql , the LPW with the same 
wave number k = 2ir/3 are created. For the stacks with 
N = 3ri + 1 (solid curves), we obtain k = 2(N - 1)-k/3N. 
This value limits to 27r/3 with an increase in N from the 
side of smaller values of k. In the stacks with N = 3n + 2 
(dash curves), we get k = 2(N + l)n/3N, which limits to 
27r/3 from the side of bigger values of k. So the idea of 
the LPW creation at the breakpoint explains the group 
behavior of IVC in Fig. 2J The value of h p depends on 
k but does not depend on N at chosen parameters a 
and /3; i.e., the creation of the same mode in the stacks 
with different N leads to the same value of Ib p . So we 
may predict a different commensurability manifestation 
in the IVC of stacks with a different number of IJJ. This 
is a generalization of the commensurability effect we have 
observed in Ref|4] at small a and (3. 

As summary, we showed that coupling between junc- 
tions changes crucially the dependence of the return cur- 
rent on a dissipation parameter. Particularly, it leads to 
the appearance of the plateau on the /3-dependence of 
the BPC on the outermost branch and the oscillation of 
the BPC as a function of (3. Using the idea that at the 
breakpoint the parametric resonance is approached and 
a longitudinal plasma wave is created, we modeled the a- 
and /3-dependence of the BPC and obtained good agree- 
ment with the results of the numerical simulation. We 
demonstrated that the study of the a- and /3-dependence 
of the BPC for the stacks with a different number of IJJ 
gives us the instrument to determine the wave number of 
the LPW. 
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